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ABSTRACT 

Improving the quality of positron emission tomography (PET) images, affected by low resolution and high 
level of noise, is a challenging task in nuclear medicine and radiotherapy. This work proposes a restoration method, 
achieved after tomographic reconstruction of the images and targeting clinical situations where raw data are often 
not accessible. Based on inverse problem methods, our contribution introduces the recently developed total gener¬ 
alized variation (TGV) norm to regularize PET image deconvolution. Moreover, we stabilize this procedure with 
additional image constraints such as positivity and photometry invariance. A criterion for updating and adjusting 
automatically the regularization parameter in case of Poisson noise is also presented. Experiments are conducted 
on both synthetic data and real patient images. 

Index Terms — PET imaging, total generalized variation, deconvolution, Poisson noise, inverse problem 

1. INTRODUCTION 

In radiotherapy, positron emission tomography (PET) images provide oncologists with useful information about 
the metabolic activity of the patient. The accumulation of the 18 F-fluorodeoxyglucose ( 18 F-FDG) radiotracer in 
tissues of abnormally high metabolic activity, e.g., cancer cells, leads to the emission of positrons later detected by 
the PET imaging system [1], Combined with anatomical CT or MR images, these functional images are widely 
used for diagnosis, monitoring during the treatment and follow-up of the patient after therapy. 

A challenging step in radiotherapy treatment is the accurate delineation of the tumour volume based on PET 
images [2], Some methods rely on visual interpretation or on activity threshold determination but suffer from 
methodological limitations [2], The use of more complex segmentation algorithms is impeded by the two main 
drawbacks of PET imaging: low spatial resolution and high level of noise [1], The first one is mainly due to the 
blur introduced by the random positron travel, the photon diffusion in the patient tissues before annihilation and 
the large size of the scintillators required to detect high-energy photons. The resulting point spread function (PSF) 
of the PET system is generally not uniform and can vary slightly in the field of view (FOV). The latter, i.e., the 
noise in raw data, is due to low photon detection efficiency of the detectors and the limitation of the injected tracer 
dose for obvious radio-protective reasons. The noise is Poissonian in the projection space, before reconstruction. 

Improving the quality of PET images is essential before applying more advanced segmentation techniques. 
Deblurring and denoising can be applied either during the reconstruction process or in a post-processing phase. 
The second approach is more appropriate for clinical use since the access to raw data or to the reconstruction 
algorithm of the scanner is not always possible. After sinogram correction and standard iterative 3D reconstruction 
algorithms like OSEM [3], noise in image is still multiplicative in first approximation (scaled Poisson). Nowadays, 
the most widely used tool in clinical denoising is a classical Gaussian filter. Edge-preserving filters such as bilateral 
filters or M-smoothers are also used but, like the Gaussian filter, they are generally designed for additive noise [4], 
This issue can be overcome using first a variance-stabilizing transform (VST) such as proposed by Anscombe [5], 
For the deblurring step, classical deconvolution algorithms like Van Cittert’s or Landweber’s [4] are used. Other 
methods combine the denoising and deblurring steps by solving an optimization problem (inverse problem) [6,7]. 

In this paper, we choose the inverse problem approach in a post-processing phase and propose a problem 
formulation specific to the restoration of PET images. Our method introduces the total generalized variation 
(TGV) as a regularization term [8] for PET deconvolution, considers specific properties of PET physics such as 
positivity and photometry invariance and takes into account the Poisson statistics of noise through the data fidelity 
term and the choice of the regularization parameter A. A criterion for automatic selection of A is presented. The 
validation of the TGV algorithm is made on synthetic images and on real medical images where the PSF is well 
estimated. 
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The paper is organized as follows. Sec. 2 introduces the mathematical model and the underlying assumptions, 
as well as theoretical and numerical aspects related to the convex optimization problem. The experimental setup is 
described in Sec. 3, along with the results, which are discussed in Sec. 4. Conclusions are presented in Sec. 5. 

2. FORWARD MODEL AND INVERSE METHOD 

In this section, we describe the forward model and the assumptions we made. We then explain the formulation of 
the inverse problem based on a maximum a posteriori estimator (MAP). Finally, we give some details about the 
numerical implementation and the choice of the parameters 

2.1. Forward model 

For the sake of simplicity, we work with a discretized model. The two dimensional images of size Ni x 
belong to the Euclidean space with N = N-\ . Let u 0 £ R w be the original and unkown PET image of the 

functional processes in a patient body. After tomographic reconstruction, the observed image z £ R w is associated 
with Uq through the forward model 

2 = Af„(K(uo)), 

where K is a blur operator accounting for both the physical and instrumental inaccuracies and A f v is a noise 
corruption of parameter v. Some assumptions on K and A f v lead to a simpler model. In first approximation, the 
PSF near the center of the FOV is uniform, Gaussian and isotropic in the 2D-plane [4], Consequently, the blurred 
image is assumed to result from the convolution of the original image with the PSF of the scanner. As mentionned 
in Sec. 1, noise is considered as Poissonian in first approximation. A simpler model is 

z = V(Ku 0 ), (1) 

where the convolution operator K £ K ,vx iV is assumed to be known (see Sec. 3), linear (uniform PSF) and 
bounded. Each pixel i of the original image is corrupted by Poisson noise V with mean (Ku (l ),. 

Physics of PET imaging suggests two additional constraints to model (1). The first property is the positivity 
since the original image u o represents a nonnegative activity in a phantom or in vivo, e.g., the metabolic activity. 
Hence, 

u 0 R 0. 

The second property is photometry invariance. Total photometry preservation means that the total photon counts 
in the original and observed images are approximately the same. This property is particularly valuable for quan¬ 
tification aspects like the measurement of standardized uptake values [1], 

N N 

i—l i—1 


2.2. Inverse problem formulation 

Solving an inverse problem consists in finding an estimator it o of the original image Uq from observations 2 in 
(1). Noise and the low-pass filter effect of the PSF lead to an ill-posed problem for which direct inversion is not 
possible and unicity of the solution is not guaranteed. To address this issue and reduce the set of possible solutions, 
a widely used method is to regularize the problem [9]. A penality term encourages the solution to respect a certain 
prior model of the original image. 


2.2.1. A priori knowledge on the PET images 

A possible assumption about the image concerns its piecewise smoothness. In many applications, frameworks 
based on total variation (TV) are used to promote piecewise constant images [8]. However, as Knoll et al. [10] 
mentionned in the case of MRI, this type of regularization is not well adapted to real medical images. Since 
these images are not piecewise constant, staircasing artifacts are observed in smooth image areas with TV regu¬ 
larization [8,10]. The total generalized variation (TGV), introduced by Bredies et al. [8] can be considered as the 
generalization of TV to higher-order image derivatives. Simplifying its presentation to our conventions, the second 
order TGV of x £ M. N is defined in a discrete setting as 

TGV^(at) = min ||Vx - w\\ 2 ,i + cv||er(tu)|| 2 ,i, (2) 

weR Nx2 
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where a € K. is a positive constant balancing between the edge-preserving term and the smoothness-promoting 
term. The general discrete gradient operator V applicable to N x k tensor fields is defined as 

V : R Nxk -> R Nx2k , x h> (Vx) = (Via:, V 2 x), 

with k € No and Vj € R NxN the first spatial derivative of the tensor field in direction e, : . The symmetrized 
derivative operator e is 


£ : R Nx2 -> R Nx 4 , x ^ e(x) = \ ((Vs) + (Vx)5 23 ), 

where S23 G {0, l} 4x4 is a matrix permuting the second and the third columns of (Vx). Applying e on the 
gradient of an image provides information about its second derivative. The concept of symmetrized tensors is 
detailed in [8], Finally, the L pq norm of x G R Nxk is defined as 

/ N 

11*11™= Em 

\*=i 

Deriving '\'GV 2 t (x) is not as easy as TV(s) because an additional minimization problem has to be solved. The 
optimal value of w in (2) depends on the features of image x. Locally, in smooth areas, w is close to Va; to give 
more importance to the smoothness-promoting term. In regions near image edges, w is close to 0 to preserve sharp 
edges like TV. These properties make TGV 2 more adapted to real PET medical images than TV. Moreover, the 
absence of staircasing artifacts improves the efficiency of subsequent segmentation algorithms. 



2.2.2. Bayesian approach 

To take into account this a priori knwoledge on Uq, the best estimator choice is the Bayesian MAP estimator, 
defined as 


Uq = arg max p(u\z) 

ue r n 

= arg min — logp(z\u) — logp(tr). (3) 

«ei K 

The first term of (3) is the fidelity term. This term encourages the estimated image to be close to observations 
z and depends directly on the noise statistics. In the case of PET images, it is the negative log-likelihood of the 
Poisson distribution, i. e ., 

N 

-\ogp(z\u) = J2(Ku-z-f(Ku)) l+ r(z), 

i—l 

where r{z) depends only on z, / is applied componentwise on vector with f(t) = logi if t > 0 and 0 otherwise 
and • denotes the Hadamard (or elementwise) product. The second term in (3) is the regularization term promoting 
the respect of the a priori knowledge, i.e., a minimal TGV, 2 , . Indicator functions n and %c take into account 
positivity and photometry invariance constraints. They equal zero if the constraint is respected and +00 otherwise. 
The nonsmooth convex optimization problem reads 

N 

min A ^{Ku - z ■ f(Ku))i + ||V« - tu|| 2 ,i + a||£(ut)|| 2 ,i +i r n(u) + i c {u), (4) 

DjeB wx2 i—l 

where A > 0 is the regularization parameter, i.e., the trade-off between the data fidelity and regularization terms. 
The value a = 2 is used in practice [10]. 


2.3. Numerical implementation 

The discrete setting is described in details in [8], Discrete operators Vi and V 2 are approximated using forward 
and backward first-order finite differences. The convolution of the image with operator K is made with a FFT. 
The undesired border effects are avoided using a “padding” method presented by Almeida et al. [11], 
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Algorithm 1 for TGV denoising and deblurring of PET images. 

1: initialization: n = 0 ; u (0) = u (0) = y £ R N ; io (0) = w (0) = 0 G R JVx2 ; p (0) = 0el*; q (0) = 0 € R JVx2 ; 
r (0) = 0 £ R wx4 ; choose r (0) = ct (0) = 0.9/||L|| 2 . 

2 : repeat 

3 : p^+P = prox a (n) F * {p (n) + a {n) Ku^) 

4: q (n+1) =prox (T („ )i ,*(q (n) + a {n) {Vu (n) - w (n) )) 

5: r (n+1 ' ) — prox CT („) F * (r*- 11 - 1 + a^e(w^)) 

6: u (n+1) = prox (n)G (M (n > + r (n) (div[q (n) ] - K*p (n) )) 

7: w^ +1) = w (n > - r (n) (-q (n) - div[r (n) ]) 

8: = 2u (n+1) - n (7l) 

9: «A n +i) = 2io (n+1) - to (n) 

10: until convergence of u 


2.3.1. Chambolle-Pockprimal-dual algorithm 

Let L : X — > y be a linear continuous operator with a norm defined by ||L|| 2 = max{||.LsE|| 2 | x £ 
X with ||cc|| 2 < 1} and G : X — > [0, +oo[ and F* : y — ► [0, +oo[ be proper, convex, and lower-semicontinuous 
functions. The Chambolle-Pock (CP) primal-dual algorithm [12] is designed to solve the following saddle-point 
problem 

min max (Lx, y) — F*(y ) + G(x), 

xGX y^y y ' 

which is a primal-dual formulation of the primal problem min F(Lx) + G(x). 

x£X 

The CP algorithm belongs to the family of proximal algorithms [13]. Such algorithms are based on the notion of 
proximal operator and can deal with optimization of non differentiable functions. Let ip be a lower semicontinuous 
convex function from X to ]— oo, +oo[ such that dom/ is non empty. The proximal operator of ip : X — > X 
evaluated in x £ X is unique and defined as [13] 


P r ox Jx) := xtgmin -\\x - x\\l + p(x). 


X(zlX 


(6) 


In our case, the presence of the two indicator functions leads to a nonsmooth and non differentiable objective 
function. The use of a proximal algorithm and particularly CP is appropriate [14]. A compact formulation of (4) 
is given by 

min Fi(Ku) + P 2 (Vm — w) + F 3 (e(w)) + G(u ), 

ue r n , 

■weR Nx2 

with straightforward definitions of Lj, F 2 , F 3 and G from (4). 

A primal-dual formulation similar to (5) of this problem is derived in a product space [15] using duality prin¬ 
ciples [9,12] and leads to Algorithm 1, i e.. the CP algorithm adapted to TGV denoising and deblurring of PET 
images. Function p* is the Legendre-Fenchel conjugate of any function p. The divergence operator div is defined 
as equal to —V* with 

V* : R^x2fc R ATxfc ; x = ( a . lj a-2) (V*®! + V 2 X 2 ). 


2.3.2. Proximal operators 

Let x £ R N . The proximal operator of the primal function G is given by a combination of positivity (max) and 
photometry invariance (average) proximal operators [13], 


prox CTG (*) 


N 


max(x, 0) - — ( m ax(*, 0)) - z) i . 


i—1 

From definition (6) and the Moreau decomposition property [13], the proximal operator of F* is, for each 
component i of x £ R N [14], 


prox^ p* (xi ) = ^ (xi + A - sj(£i - A ) 2 + TtXtT) , 

if z, is non-zero and is equal to parameter A otherwise. 

By conjugation of F> and F 3 , F 2 and F 3 are indicator functions. Their proximal operators are projection 
operators onto the convex sets Q = {x £ M JVx2 |||:r|| 2j00 < 1} and R = {x £ K Arx4 |||a:|| 2i00 < q-}. 
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Fig. 1 : Results on synthetic data. From left to right: original image, corrupted images (convolution with isotropic Gaussian 
kernel and Poisson noise: j3 lop = 1 and /^bottom = 0.1), TV restored image (A t0 p = 32.8, Abottom = 4.7), TGV restored image 
(Atop = 25.3, Abottom = 5.6) and residual image computed as the difference between the original and the TGV restored images. 


2.3.3. Parameters choice 

Step-sizes er and r are updated at each iteration depending on the size of the primal and dual residuals [15]. To 
ensure the convergence, condition <jt||L ||2 < 1 has to be verified [12]. 

The choice of the regularization parameter is based on the discrepancy principle for Poisson noise adapted 
in [16] to images with null background. Parameter A is selected such that 

KLy, (7) 

where M < N is the number of non zero pixels and KL is the Kullback-Leibler divergence defined as 

N 

KL(at, y) = y -x + x- f(x) - x ■ f{y))i. 

i -1 

Interestingly, KL(z , Ku\) = F\(Ku) — Fi(z), so that (7) can be interpreted as a distance in the range of F\oK. 


3. EXPERIMENTS 

Synthetic images. The synthetic image used for a first validation of the method is a modified Shepp-Logan 
phantom with intensities in [0, 255] (see Fig. 1). The original image is piecewise constant and is therefore hardly 
representative of real medical data. In the modified version, constant surfaces are replaced with a mixture of affine, 
Gaussian and sinusoidal functions. Since the background is expected to be null in PET, it is set to zero in the 
phantom image. 

To simulate different acquisition times, the pixel intensities of the original phantom Uq are multiplied by a 
constant positive factor (3 {13 values in [10 —2 ,10 2 ]). The bigger j3, the higher the number of photon count per 
pixel and the smaller the relative effect of the noise. To simulate the blur, the 128 x 128 images are convolved with 
a normalized Gaussian kernel with standard deviation equal to 1.17 pixel. This PSF is assumed to be ideal, i.e., to 
be constant over the whole FOV. Finally, each pixel of the image is corrupted by Poisson noise. Algorithm 1 was 
used to deblur and denoise the images with TGV regularization (as well as TV with modification of the proximal 
operators). Parameter A is updated at each of the 20 meta-iterations according to the following rule: 

x(I+l) _ 

M/2 ’ 

where the updating factor is called “KL ratio”. The signal-to-noise ratio (SNR) is used to quantify the quality of 
the corrupted image (SNRj„) and the restored image (SNR 0 , (i ) relative to the original one. Convolution and noise 
influence the SNR,,,. SNR between image x and the original image y is defined as 

S NR(x, s) = 201„g lo (Jhh E ). 

Results. A visual comparison between TGV and TV regularizations is presented on Fig. 1 for two levels 
of noise: (3 = 1 and (3 = 0.1. The impact of the noise level (or similarly, the acquisition time) on the SNR 


5 




• TGV 
oTV 


20 


SNR out [dB] 



3.4 dB g 



10 


10.4 dB 


-2 



4> log 10 p [-] 
2 


Fig. 2: Output SNR as a function of the level of noise @ for TGV and TV regularizations. Input SNR are indicated for extremal 
values. 


is illustrated in Fig. 2 for both TGV and TV regularizations of the same corrupted images. The input SNR is 
indicated for extremal values. Finally, Fig. 3 shows the evolution of the KL ratio as a function of the meta-iterations 
for automatic selection of parameter A. 

Patient images. PET images of patients with pharyngolaryngeal squamous cell carcinoma were acquired on 
a Siemens Ecat HR scanner for a previous study [17]. The size of the image is 128 x 128 x 47 and the voxel 
size is 2.2 x 2.2 x 3.12 mm 3 . The algorithm was applied on axial slices. The PSF was measured experimentally 
with a line source perpendicular to the axial slices in diffusing material (water). In first approximation, the PSF 
at 100 mm from the FOV center was Gaussian, isotropic, with 6 mm of FWHM. Algorithm 1 was used to deblur 
and denoise the images. The value of the optimal regularization parameter was set to A = 5 (see Sec. 4). 

Results. Since there is no original image available, main results consist of a visual comparison between TGV 
and TV regularizations (see Fig. 4). 


4. DISCUSSION 

Synthetic images. Results on Fig. 1 show that both methods remove the noise and preserve edges. Some stair¬ 
casing artifacts affect TV restored images since the original image is not piecewise constant. This cartoon-like 
appearance removes details like the small white disks. As expected, TGV restored images keep the smoothness 
property of the original image. The value of parameter A depends on the image dynamics (bigger if higher in¬ 
tensity) and is similar for the two regularization methods. Distinction between TGV and TV is not so obvious 
on the graph of the SNR evolution. In average, SNR out of TGV is 0.5 dB above TV’s and is mostly due to the 
contribution of small (3 (2 dB for /3 = 10 -2 ). The updating rule for parameter A gives convincing results since 
the KL ratio converges to 1 and the standard deviation to 0 (average on 13 realizations). These results validate 
the automatic selection rule of A in case of Poisson noise. Regarding Fig. 3, 10 iterations seem to be enough for 
accurate estimation of A. 

Patient images. The staircasing artifacts observed in synthetic image with TV restoration are also visible on 
real medical images (Fig. 4). Since data are not corrupted with a pure Poisson noise (see Sec. 1), the updating rule 
for regularization parameter A cannot be applied. The choice of the value of A is highly subjective and depends on 
the application as mentionned in [14]. 

5. CONCLUSIONS AND FUTURE CHALLENGES 

This work presents a new method for deconvolution of PET images based on TGV regularization. The use of 
TGV instead of TV as a prior is more adapted to real medical images for which the piecewise-constant assumption 
is not verified. In case of pure Poisson noise, we validate an updating rule for parameter A. Further work will 



Fig. 3: KL ratio evolution as a function of the meta-iteration count. The thick black curve indicates the mean of 13 realizations 
(different noise levels) and the thin curve indicates the standard deviation. 
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Fig. 4: Results on real medical data (full and zoom). From left to right: corrupted PET image, TV and TGV restorations with 

A = 5. 

investigate the nature of the noise after reconstruction (with Monte-Carlo simulations) to improve the formulation 

of the inverse problem as well as the updating rule of A. 
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